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ABSTRACT 

The 13Myr old star HD 106906 is orbited by a debris disk of at least 0.067 Mmooh 
with an inner and outer radius of 20 AU and 120 AU, respectively, and by a planet at 
a distance of 650 AU. We use this curious combination of a close low-mass disk and 
a wide planet to motivate our simulations of this system. We study the parameter 
space of the initial conditions to quantify the mass loss from the debris disk and its 
lifetime under the influence of the planet. We find that when the planet orbits closer 
to the star than about 50 AU and with low inclination relative to the disk (less than 
about 10°), more disk material is perturbed outside than inside the region constrained 
by observations on timescales shorter than 1 Myr. Considering the age of the system, 
such a short lifetime of the disk is incompatible with the timescale for planet-planet 
scattering which is one of the scenarios suggested to explain the wide separation of 
the planet. For some configurations when the planet orbit is inclined with respect to 
the disk, the latter will start to wobble. We argue that this wobbling is caused by a 
mechanism similar to the Kozai-Lidov oscillations. We also observe various resonant 
structures (such as rings and spiral arms) induced in the disk by the planet. 

Key words: celestial mechanics - planet-disc interaction - planetary systems: for¬ 
mation - circumstellar matter - planets and satellites: individual: HD 106906 b - open 
clusters and associations: individual: Lower Centaurus Crux 


1 INTRODUCTION 

About a dozen planetary mass companions at wide sep¬ 
arations of about 50-100 AU from their host stars have 
been revealed by direct imaging surveys during the past 
decade (Kraus et al.||20i4| and several cases were observed 


at separations of 150-300 AU (e.g., Lafreniere et al. 2008 
[Kraus et al.|[2014[ ). Moreover, two recent discoveries re¬ 
port companions located as far as ~650AU (Bailey et al. 


2014) and ~2000 AU ( Naud et al.|2014 1. The origins of such 


wide planetary mass companions is not well understood and 
presents important constraints for our general understand¬ 
ing of planet formation. Several scenarios have been pro¬ 
posed, and depending on the eccentricity and separation of 
the planet, environment in which the system evolves, and 
timescales of the formation, two main mechanisms are usu¬ 
ally considered. 


In situ formation by core accretion (e.g., Ra fikov|2011 
or protoplanetary disk fragmentation (e.g. 


Boss 2011 


Vorobyov 20131 can explain part of the observed popula¬ 


tion of the wide orbit planets but is unlikely to be the only 
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formation channel (see also Veras et al. 2009 D’Angelo et al. 


2011 or Chabrier et al. 2014 for recent reviews of the topic). 


Another explanation argues that the planet formed 
closer to the parent star in the protoplanetary disk and 
was scattered outward by dynamical interaction with an¬ 
other planet system or with perturbation of external origin 
(see e.g., Davies et al. 2013| for a summary on various in¬ 
teractions in planetary systems). Given the diversity of the 
observed wide planetary systems and the environment they 
are expected to form in, the parameter space for the ini¬ 
tial conditions of such scattering events is extremely large 
and complex. The formation can involve for example, stellar 
flybys (e.g., Malmberg et al. 2011), exchange interactions 
( Portegies Zwart fc McMillan||2005| ), planetary migration 
(e.g., Crida et al.||2009 ) and scattering in a multiple plane¬ 
tary system ( Scharf fc Menou||2009 1, dynamical interaction 
between circumstellar disks and planets (seejBaruteau et al.j 


2013 


for a recent summary), the effects of Galactic tides 
(e.g., Veras fc Evans||2013|), recapture of free floating plan¬ 
ets (jPerets &; Kouwenhoven||2012|), or combination of these 


interactions (Raymond et al. 2010 Boley et al. 2012 Hao 


et al. 2013). Studying specific objects narrows down this 


parameter space since some of the characteristics are con¬ 
strained by observations. 
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In this context, we focused on HD 106906 which is a 


F5-V star with a debris disk (Chen et al. 2005 20111 and 
a planetary mass companion at a distance of about 650 AU 
I Bailey et al.|2014 1. The chance of coincidental projection of 
the star and planet is negligible, and therefore the observed 
distance between the star and the planet is interpreted as 
an orbital separation. Irrespective of the inclination of the 
planetary orbit, which is unknown, the observed separation 
must be part of the orbit, which makes it one of the widest 
separation ever observed. 

Regardless of the process that caused this planet to have 
such a wide orbit, the observed debris disk has survived. The 
lifetime of the debris disk as is observed, constrains how long 
ago the current conhguration formed. In this paper we study 
the timescale on which the disk erodes due to the influence 
of the planet, and use this timescale to constrain the mech¬ 
anism that delivered the planet in its extremely wide orbit. 
We carry out simulations of the evolution of the disk under 
the influence of the planet, taking the observed character¬ 
istics of the system as the initial conditions. We vary the 
inclination of the disk with respect to the planetary orbit 
and the pericenter distance of the planet (i.e., its eccentric¬ 
ity under the assumption that the apocenter distance of the 
orbit is 650 AU) within the observational constraints, and we 
explore the erosion timescale of the disk due to the planet. 


Table 1. Characteristics of the HD 106906 system. 


Characteristic 

Value 

Unit 

Ref. 

Distance 

92±2 

pc 

a 

Age 

13±2 

Myr 

b 


HD 106906 A 



Spectral type 

F5V 


b 

Mass M* 

1.5±0.1 

M© 

b 

Luminosity L* 

5.6±0.8 

L© 

b 

Temperature 

6516±165 

K 

b 


HD 106906 b 



Mass Mb 

11±2 

Mjup 

c 

Separation Hb 

650±40 

AU 

c 


disk 



24 pm flux density 

103.1±2.5 

mjy 

d 

70 pm flux density 

281±28 

mjy 

d 

Fractional luminosity Ljr/L* 1.3 x 10 ^ 


d 

Dust grain temperature 

95 

K 

c 

Inner radius 

~20 

AU 

c 

Outer radius 

<120 

AU 

c 

Minimum mass 

0.067 

Mjvioon 

d 


References: a -|van Leeuwen] ||2007|l , 6-|Pecaut et al.| ( |2012^ ; c - 
[Bailey et al.| l |2014| |; d- jChen et al.j | |2011| |. 


1.1 The HD 106906 system 

HD 106906 (or also HIP 59960) belongs to the Lower Centau- 
rus Crux (LCC) group which is a subgroup of the Scorpius- 
Centaurus (ScoCen) OB association | |de Zeeuw et al.[jl99'^ . 
The host star, called HD 106906 A, is classihed as F5-V star. 
Pecaut et al. ( |2012| ) measured the median age of the LCC 
group of 17 ± 1 Myr, and the mass and age for HD 106906 A 
of M* = 1.5 Mq and 13 ± 2 Myr, respectively. In Table ^ 
we summarize the observed data and derived characteristics 
of the HD 106906 system. 

The observed infrared (IR) spectral energy distribution 
of HD 106906 A shows a strong excess that indicates the 
presence of a debris disk with inner cavity.|Chen et al.|||2011[ 
see also [Chen et al.|2005| for the initial results based on the 
same observational data) obtained broadband observations 
of HD 106906 with the Multiband Imaging Photometer for 
Spitzer at 24 and 70 pm. By htting these excess fluxes with 
a single black-body, they derived the disk’s color temper¬ 
ature of 93 K and fractional IR luminosity with respect of 
the star Lir/L* = 1.3 x 10“^. Bailey et al. (20141 confirmed 
these results using additional Spitzer data up to ~100 pm, 
obtaining a disk temperature of 95 K. 

The disk around HD 106906 A is expected to be opti¬ 


cally thin. Chen et al. (2011 \ identified 55 stars with IR ex¬ 


cess in their sample of 167 ScoCen OB Association members 
of intermediate-age (10-30 Myr) and F-, G-, or K- spectral 
types. They did not find any significant difference between 
the distribution of the IR excess (measured by the Lir/L* 
ratio) for fast and slow rotating stars. As a difference is 
expected in rotation speed for stars hosting gas-rich and 
gas-poor stars (due to magnetic braking, e.g., [Rebull et ah] 
2006), it is likely that the stars in ScoCen association have 


optically thin and gas poor disks. 

Since the disk is not resolved at any wavelength, its 
characteristic extent can be estimated from the tempera¬ 


ture. Assuming the dust grains are black-bodies in radiative 
equilibrium with the central star, an optically thin disk with 
grains of constant size and chemical composition, [Chen et al.j 
(20111 derived a single grain distance of about 34 AU. Based 


on the comparison with Hershel observations of a sample of 
resolved circumstellar disks, [Bailey et al.j ( |2014| ) further es¬ 
timated the extent of the disk to be about 20-120AU (for 
the optically thin disk). Chen et al. (2011) also estimated 
the minimum dust grain size of 1.4 pm, and the minimum 
mass of the IR-emitting dust grains of 0.067 Mmooh- 

The planetary mass companion of HD 106906, called 


HD 106906 b, was discovered by Bailey et al. (20141 with the 


Magellan Adaptive Optics/Clio2 system. They obtained re¬ 
solved images of the companion, confirming that the planet 
is comoving with the host star, and classified its spectral 
type as L2.5±l. As mentioned above, the projected sep¬ 
aration between the host star and the companion then is 
650 AU. Using evolutionary models for an object of this spec¬ 
tral type and age corresponding to the one of the LCC group, 


Bailey et al. (20141 further estimated the mass of the planet 


to be Mb = ll±2Mjup. Properties of the planet make the 
formation of HD 106906 difficult to explain. The two most 
compelling formation mechanisms for the origin of planets 


in wide orbits are discussed by Bailey et al. (2014): i) in situ 


formation at a large separation, as wide as the orbital sep¬ 
aration found in some binary stars; and ii) formation in a 
tight orbit and the subsequent scattering to the current wide 
orbit. The mass ratio Mb/M* ~ 0.01 is unusually small for 
the first suggested mechanism. In the later scenario a per- 
turber must have been present in order to move the planet 
to its current orbit. This culprit however, may be long gone, 
lost in interstellar space. This is consistent with the lack of 
another massive planet in the system ( Bailey et al.|20f4 l — 
no other object is detected within the observational limits 
which translate to a mass no grater Mb beyond 35 AU, and 
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a mass no greater than 5-7Mjup beyond 70 AU. We cannot 
rule out other formation mechanisms, such as the possibly 
capturing of the planet from the surrounding environment 
in the LCC group. 

Here we explore the lifetime of the current configura¬ 
tion of the system. Planet-planet scattering is predicted to 
occur after the dissipation of t he gas from the circum stel- 
lar disk at about 10® yr (see e.g. Chatterjee et al.|2008 and 


an extension of the mixed variable symplectic scheme, which 


references therein). Planets at wide separation (> 100 AU) 
are estimated to be most probably produced on timescale s 
up to 10^ yr (e.g., Veras et al.|2009 Scharf fc Menou|2009 |. 
If the current planetary orbit is the result of a scattering 
interaction with another planet, both planets once orbited 
the parent star in a much closer orbits, probably within the 
observed inner edge of the disk. The current planetary orbit 
must still bear the memory of that original orbit and the 
place where the scattering happened, closer to the parent 
star, should also be part of the orbit. The lifetime of the 
disk under the influence of such a planet should then be at 
least a few Myr in order to be consistent with the lifetime 
of the system. 

We investigate the mass loss of the disk for different 
eccentricities and inclinations of the orbit with respect to 
the disk. 


2 SIMULATIONS 

We performed simulations of the evolution of the system 
starting with initial conditions corresponding to its current 
observed characteristics (see Table [^. We varied some of 
the unconstrained properties, namely the pericenter of the 
planetary orbit, Rp, and the inclination of the disk, i, since 
these can in principle be random depending on the formation 
process of the system. 


2.1 Method 

We calculated the orbit of the star-planet system inde¬ 
pendently of the evolution of the disk. Since the mass of 
the disk is small compared to the planet or the star, we 
represented the disk by a number of zero-mass particles — 
planetesimals — and hence we do not take the self-gravity of 
the disk into account. 

All calculations were carried out within the Astrophysi- 


cal Multipurpose Software Environment or AMUSE (Porte- 


gies Zwart et al. 

2009 

Pelupessy et al. 

2013f 


A-body integrator HUAYNO ( [Pelupessy et al.|2012[ ) to cal¬ 
culate the orbit of the star-planet system. The orbits of the 
disk particles were calculated by solving the Kepler’s equa¬ 
tions using universal variables (adopted from the SAKURA 
code, Gongalves Ferrari et al.|2014 (. The implementation of 
the solver in AMUSE allows us to efficiently integrate Kep- 
lerian orbits in the potential of a central star with a number 
orbiters (i.e., planetesimals orbiting the star in our case). 
Our approach is not self consistent — the planet and the 
star are not influenced by the planetesimals in the disk. The 
gravitational influence of the planet is coupled with the plan¬ 
etesimals. This coupling, called Bridge (Fuji! et al.|2007|, is 


was developed by Wisdom & Holman (19911, and it is used 


here to couple different dynamical regimes within one self- 
gravitating system (i.e., the planetesimal debris disk and the 
planet orbiting the central star). The time complexity of our 
numerical scheme is oc A, rather than the usual oc A^ for a 
direct A-body approach. The implementation of Bridge in 
AMUSE is described in |Pelupessy et ah (|2013 (. 

The symplectic mapping method of Wisdom fc Holman| 


119911 was first applied to calculate the long-term evolution 


the solar system and has since been widely used to simulate 
the evolution of planetary systems in general, including in¬ 
teraction with planetesimals. Fragmenting planetesimals are 
generally considered to be the parent bodies of the dust that 
is observed as a debris disk (e.g., Wyatt|2008 1 and complex 
methods have been developed to accurately model this pro¬ 
cess (see, e.g., Thebaui^|2012 and references therein). The 
planetesimal disk approximation is often used to define the 
spatial and velocity distributions of the dust particles. For 


example, Larwood & Kalas (20011 investigated the affect of 


stellar flybys on the structure of the debris disk observed in 


the P Pictoris system, and similarly in Chiang et al. (20091 


for the Fomalhaut system. Wyatt (20031 or Reche et al. 


(20081 studied the resonant trapping of planetesimals due 


to planetary migration. Lestrade et al. (20111 investigated 


the stripping of the planetesimal debris disk by a close stel¬ 
lar flyby. Long-lived asymmetric structures were simulated 
by, e.g., Faramaz et al. _^2013|eccentric debris disk around 
Cf Reticuli) or Pearce & Wyatt (2014 more general case of 
a planet within the outer edge of the disk). 

We tested the method by comparing our implementa¬ 
tion with direct A-body integrations, which gave qualita¬ 
tively and quantitatively the same results; and we success¬ 


fully reproduced the results of Lestrade et al. (20111. 


http://amusecode.org 


2.2 Numerical setup and initial conditions 

Following the observations, we assumed a mass of 1.5 Mq 
for the star and llMjup for the planet (see Sect. o and 
Table [^. The apocenter distance of the planet was 650 AU 
in all our simulations. This is the observed separation, which 
we assume to be the apocenter of the orbit, and which is the 
planet’s initial position in our simulations. The pericenter 
distance of the planet, Rp, had values ranging from 1 AU to 
650 AU, corresponding to orbital eccentricities of 0.997 and 
to circular orbit, respectively (see Table for the list of all 
pericenter values considered). The orbit of the planet was 
integrated with HUAYNO using the HOLD drift-kick-drift 
integrator. The HUAYNO integrator uses individual time- 
steps that are proportional to inter-particle free-fall times 
and the coefficient of the proportionality is called rj. We 
chose different values of rj for different pericenters (i.e., or¬ 
bital eccentricities) so that the energy conservation of the 
star-planet system is always at 10“® level and lower; this 
level of energy conservation turns out to be very conserva¬ 
tive ( Portegies Zwart fc Boekholt|2014 1. The values of rj are 
specified in Table |2| for each orbital configuration. 

Disk planetesimals begin in an initially a uniform ran¬ 
dom distribution in radius between the inner and outer disk 
radii of 20 AU and 120 AU, respectively, which corresponds 
to the values estimated from observations (see Sect. o and 
Table [^. Such choice of radial distribution corresponds to 
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Table 2. Planetary pericenters and time-steps for the inte¬ 
grations. 


Rp [AU] 

V 

iBR0 

1 

0.001 

0.001 

10 

0.001 

0.002 

20, 30, 40, 50, 60 

0.001 

0.01 

70, 80, 90, 100, 110 

0.001 

0.05 

120, 150, 200, 350, 500, 650 

0.003 

0.05 


“ The Bridge time-step, teRi is given in the units of the period 
of the circular orbit at 20 AU from the star, which is 73 yr. 


the surface density profile oc 1/r, where r is the radial 
distance to the star, which is often used to model proto¬ 
planetary disks (see e.g., Steinhausen et al. 2012 and ref¬ 


erences therein) and is supported by observations (e.g., An¬ 


drews & Williams 


20071. Following the discussion in Stein¬ 


hausen et al. (20121, we tested how our results depend on the 


chosen initial surface density profile. Since the disk is rep¬ 
resented by test particles (i.e., its self-gravity is not taken 
into account), different surface density profiles can be taken 
into account in the post-processing of the simulations. We 
considered a surface density profile oc 1/r^'®, corresponding 
to the Minimum Mass Solar Nebula ( Hayashi|1981 1, and we 
found that such profile changes the disk fractions presented 
in Sec. |3.1| by less than 10%. 

The planetesimals are initially placed in one plane with 
a random, uniform azimuthal distribution and they have cir¬ 
cular orbits. The inclination of the disk with respect to the 
planetary orbit, i, has values between 0° and 180°, where 
f = 0° corresponds to coplanar prograde case, and i = 180° 
corresponds to coplanar retrograde case. The disk plane is 
rotated around axis perpendicular to semi-major axis of the 
planetary orbit. Each simulation was carried out with 10'* 
particles, but we confirmed that increasing this number to 
10® does not change the results. Decreasing the number 
to 10® particles gives qualitatively similar results, but the 
smaller number of particles makes post processing analysis 
harder due to the lower statistics. 

The planetesimals feel the gravitational force from the 
planet with specific time-step of the interaction, Ibr — the 
Bridge time-step — which is the time step in which the sys¬ 
tem integrates the combined solver. The time-step differs for 
different initial conditions of the planetary orbit — for more 
eccentric planetary orbits we adopted a shorter time-step. 
Ibr has values ranging from 10“® to 5 x 10~® of an orbital 
period of the initial inner disk edge of 73 yr (which is the case 
for the adopted 20 AU). The values of Ibr are specified in 
Table|^for each orbital configuration. We verified the choice 
of tBR by comparing the integrations using Bridge to the cal¬ 
culations where the whole system was treated by the A^-body 
code. These control A-body simulations were carried with 
10® zero-mass particles in the disk. We used the HUYANO 
integrator in AMUSE with choice of rj giving the energy con¬ 
servation of order 10“® or lower. To treat close encounters of 
the planetesimals with the star, we use Plummer softening 
with smoothing length e = 0.001 AU = 0.2 Rq. The results 
of the direct and our Bridged direct-Kepler solver are in a 
good agreement. More quantitatively, we compared the disk 


fraction , /j/ t,—the main output of our simulations defined 
in Sect. |3.l] —which generally agrees on a ~5% level. 


3 RESULTS 

In Fig. we show an example of our simulation — the con¬ 
figuration with pericenter of 20 AU (at the inner edge of the 
disk) and the disk inclination of 5°. The surface number 
density of planetesimals in the plane of the planetary orbit 
(xy) and the edge-on plane {xz) is plotted in the upper and 
lower panel respectively. As the planet plunges through the 
disk, it perturbs the planetesimals’ orbits and the disk is 
disrupted. Some planetesimals move outside the initial disk 
region and some become unbound from the star and escape 
from the system. The majority of the particles that are mov¬ 
ing outside the initial disk region are perturbed farther away 
from the star, i.e. their semi-major axis is larger than the 
outer disk radius of 120 AU (indicated by the gray ellipse 
in Fig. [^, and only a small fraction of particles are orbit¬ 
ing within the inner disk edge (with semi-major axis smaller 
then 20 AU). Note that we do not consider collisions between 
the planetesimals themselves neither with the star nor the 
planet and no particles are removed from the simulation. 


3.1 Parameter space study 

We explored the parameter space of the pericenter of the 
planetary orbit (Rp) and the inclination of the disk with re¬ 
spect to the orbital plane (i). In Fig.j^we show the fraction 
of the disk particles that stay bound to the star after 1 Myr 
of the evolution — ribound/ntot, where ribound is the num¬ 
ber of bound particles and ntot is the total number (i.e., 
ntot = 10^). Fig. I maps the prograde cases (0° < i ^ 90°); 
the results for the retrograde configurations are generally 
similar (see below for some examples). We see that only in 
the coplanar case when the pericenter is smaller that the 
outer disk radius, a substantial number of particles is lost 
(unbound) from the system. It is hardly surprising that the 
highest number of unbound particles is produced in such 
configurations, but it is interesting that more than ~ 80% 
of the particles stays bound for all the other considered con- 
hgurations during the hrst 1 Myr. 

The number of bound particles measures what part of 
the original disk is kept within the system which, however, 
does not directly correspond to the observed disk. For ex¬ 
ample, in the second and the last snapshots of Fig.[^ we see 
that a substantial number of the planetesimals is located 
outside the disk area as it was constrained from the obser¬ 
vations. Most of these planetesimals are however still bound 
to the star and the ratio ribound/ntot is about 0.8 at IMyr 
(see Fig. [^. Majority of these bound particles perturbed 
from the disk extent have semi-major axis larger than the 
outer edge of the disk of 120 AU, while only small fractions 
orbits within the inner edge. 

To estimate how consistent our simulations are with 
the observed disk, we follow the ratio of the number of par¬ 
ticles with their instantaneous distance from the star within 
the observationally constrained disk extent and the number 
of particles bound to the star. We call this quantity disk 
fraction /d/b and it is given as /d/b = n(20AU < R < 
120 AU)/nbound, where n{R) is the number of particles at 
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X [AU] X [AU] X [AU] X [AU] 


Figure 1. Snapshots from the simulation with ijp = 20 AU and i = 5°. Time of the snapshots is indicated above each panel in Myr and 
in Pu (orbital period of the planet). The color-scale maps the number of planetesimals, rip, projected in the planetary orbit plane {xy) 
and the edge-on view of the initial disk (plane perpendicular to the planetary orbit, xz) in the upper and the lower panel respectively. 
The star, the planet and its orbit are indicated in blue. The gray ellipse and line segment show the initial extend of the disk. The planet 
and the planetesimals rotate in the same sense, counter clock-wise in the xy plane. 



Pp [AU] 


Figure 2. The fraction of particles that stay bound to the star after 1 Myr mapped in the pericenter-inclination plane. The planetary 
pericenters and disk inclinations are changing along the horizontal and the vertical axis, respectively. The plane is divided in colored bins 
and the Rp and i of the used grid are indicated by points. Note that the horizontal axis is logarithmic except for the smallest pericenter 
(1 AU), which is shown in different scale for clarity. 


given distance R (spherical radius) from the star. We use 
the instantaneous distance because the disk is not resolved 
in the observations and its extent is estimated from the tem¬ 
perature that is given by the distance of the debris from 
the star. We tested that in case when the semi-major axis 
of the particles’ orbits is used instead of the instantaneous 
distance, the evolution of the ratio stays generally similar 
however, its modulations, both the short- and the long-term 
(see Sec. 1^, are not present. 

As mentioned, the ratio measures the similarity 
of the simulated system to the observed state. If this ratio 
is high, most of the particles are orbiting within the radii 
constrained by observations; low value of /d/b indicates that 
most of the particles bound to the star are orbiting outside 
the constrained radii. 

In Fig. we show the evolution of /d/b over IMyr for 


the cases when the pericenter of the planetary orbit is 1 AU 
and when it coincides with the inner edge of the disk (7?p = 
20 AU) for a number of disk inclinations. We focus on the 
cases with the pericenter within the inner disk edge because 
such configurations are expected if the planetary orbit is the 
result of a planet-planet scattering. In both cases, generally 
the lower the inclination, the lower the ratio /d/b and there 
is about 30% difference between the inclination of 5° and 
the coplanar configuration. The evolution of /d/b (I) is not 
monotonic and is subject of (at least) two modulations with 
different timescales of about 0.05 and 0.3 Myr. 

In Fig. we show /d/b (I) for configurations when the 
disk has a retrograde rotation with respect to the orbit of 
the planet (i.e., i ^ 90°) with pericenter of 20 AU. The evo¬ 
lution of the disk fraction looks generally very similar to the 
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Figure 3. Evolution of the fraction fd/hW for a pericenter distance of 1 AU (left) and 20 AU (right) and various inclinations of the disk 
with respect to the planetary orbit i < 90° (prograde cases). The lines of different colors correspond to different i as indicated to the 
right of each plot. 



t [Myr] t [Myr] 

Figure 5. Evolution of the fraction /d/b(^) for disk’s inclination of 0° (left) and 45° (right) and various pericenters of the planetary 
orbit. The lines of different colors correspond to different Rp as indicated to the right of each plot. 



t [Myr] 


Figure 4. Evolution of the fraction /d/b(^) for a pericenter dis¬ 
tance of 20 AU and various inclinations of the disk with respect 
to the planetary orbit i ^ 90° (retrograde cases). 


prograde cases with the same planetary pericenter (Fig. 
right). 

Finally, in Fig. we show /d/b(^) for fixed inclinations 
of 0° and 45° and several values of the pericenter of the 
planetary orbit. As expected, the disk fraction is generally 
higher for the confignrations with larger pericenters — more 
than about 80% of the particles is within the disk for pari- 


centers beyond the onter edge, Rp > 120 AU. Similarly as 
in Figs. and the disk fraction oscillates with two differ¬ 
ent timescales — the modulation with the longer timescale 
occurs only in cases with non-zero inclination, while the 
shorter one is present for configurations with higher disk 
fraction /d/b ~0.7. The possible explanation of these is dis¬ 
cussed in Sec. ID 


3.2 Disk lifetime 

When the ratio /d/b decreases below 0.5, more bound 
disk particles are located outside than inside the distance 
range constrained from observations. The moment when 
/d/b(to. 5 ) = 0.5 can be taken as a measure of the lifetime of 
the disk as we observe it today. In Fig. [^we show how to .5 
changes with pericenter Rp for different inclinations. Note 
that for some of the simulations to obtain to. 5 for pericenters 
Rp = 1 and 10 AU, 10® particles were used rather than stan¬ 
dard 10^. We tested that this does not change the results (see 
also Sec. 2.2[ ). In some configurations, /d/b(I) is not mono¬ 
tonic and the moment when /d/b = 0.5 occurs more than 
once (see Sec. for discussion on the oscillations and wob¬ 
bles) and we use the earliest moment to measure to .5 in these 
cases. Using the later times leads to qualitatively similar plot 
and does not change the conclusions. Fig. shows the to .5 
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Figure 7. The ratio /d/b S't the time of 1 Myr mapped in the pericenter-inclination plane. The planetary pericenters and disk inclinations 
are changing along the horizontal and the vertical axis, respectively. The plane is divided in colored bins and the Rp and i of the used 
grid are indicated by points. The color maps the /d/b(lAlyr) for given configuration of Rp and i. The horizontal axis is logarithmic 
except for the smallest pericenter (lAU), which is shown in a different scale for clarity. Contour lines are over-plotted, their levels go 
from 0.5 and are increasing by 0.1; the contour for /d/b(l Alyr) = 0.5 is indicated by the dashed line. 



Rp [AU] 

Figure 6. Dependence of to. 5 , when fd/hdo.s) = 0.5, on the 
pericenter of the planetary orbit for different inclinations. The 
dashed horizontal line indicates the lifetime of the system, 13 Myr. 


for pericenters up to 150 AU; wider pericenters, regardless 
the inclination, have to .5 longer than the system lifetime. 

The timescale to ,5 is shorter than 1 Myr for the config¬ 
urations with low inclination {i ^ 10 °) and the pericenters 
smaller and close to the inner edge of the disk {Rp ^ 60 AU). 

The choice of /d/b = 0.5 as the critical value to test 
for consistency with the observations is arbitrary. The ap¬ 
propriate choice is in principle given by the observational 
limits (i.e., the minimal detectable mass-density of the de¬ 
bris disk). We verified that the general results do not change 
when considering a /d/b of 0.3-0. 8 . As expected, the lower 
the ratio (i.e., the smaller the fraction of the particles within 
the original disk region) the longer the timescale. 

Values of /d/b at 1 Myr are shown in Fig.j^ Similarly as 
in Fig. more than half of the bound particles are located 
outside the disk (i.e., /d/b(lMyr) < 0.5) for the small peri¬ 
centers and the low inclinations. The disk stays relatively 
unperturbed for Rp ^ 150 AU regardless the inclination. 


4 DISCUSSION 


4.1 Disk wobbling and Kozai—Lidov-like 
oscillations 


As mentioned in Sec. |3.2[ for some of the configurations 
with inclined disks, the disk fraction does not decrease 
monotonously (see Fig. The modulation in /d/b(f) can 
be explained by a wobbling of the disk. We argue that this 
wobbling is caused by a mechanism similar to Kozai-Lidov 
oscillations ( Kozaip962 Lidov 19621. 

The Kozai-Lidov mechanism describes exchange of an¬ 
gular momentum in stable hierarchical three-body systems. 
The inner binary is periodically excited to high eccentric¬ 
ity and inclinations with respect to the initial orbital plane, 
and its argument of periapse librates (i.e., oscillates around 
a fixed value) with the same period. However, the energy, 
i.e., the semi-major axis of the orbit, does not change in 
the standard picture of the Kozai-Lidov mechanism (e.g., 
Mardling & Aarseth 20011. The amplitude of the oscilla¬ 


tions depends on the relative inclination of the orbits — the 
higher the inclination the bigger the changes of eccentricity 


(e.g., Innanen et al. 19971. The period of the Kozai-Lidov 


oscillations depends on the masses of the bodies, the periods 
of the orbits, and the eccentricity of the outer binary. 

The Kozai-Lidov timescale for the restricted three-body 
problem is approximately given by (see, e.g., [Hamers et S(] 
2013 and references therein). 


Tkl = a 


-b Mb 

Pd Mh 




2 \ 3/2 


( 1 ) 


where Pb and es are the period and eccentricity of the plan¬ 
etary orbit, respectively. M* and Mb are the central star 
and the planet mass, respectively. The orbital period of the 
disk planetesimal is Pd. a is a coefficient of order unity. 

The strongest modulation of fd/h{t) in Fig. happens 
for the case with Rp — 1 AU and i = 45°. This configuration 
(nor the others presented in Figs.j^ does not correspond to 
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Figure 8. Snapshots from the simulation with Rp = 1 AU and i = 45°; see Fig. for detailed description. 


the classical Kozai-Lidov example — the planet orbits inside 
the inner disk radius and the system star-planet-disk par¬ 
ticle does not classify as hierarchical triple. However, since 
the planetary orbit is very eccentric (eccentricity of 0.997 
for Rp = lAU), the time the planet spends closer to the 
star then 20 AU, is extremely short — less than 0.3% of the 
orbital period — and the time within the outer disk radius of 
120 AU is about 3.6% of the period. The planet moves out¬ 
side the disk for most of the time and periodically perturbs 
the orbits of the disk particles, changing their inclination 
and eccentricity similarly to the Kozai-Lidov mechanism. 
At the same time, we do not observe substantial change in 
the semi-major axes of planetesimals’ orbits and the modu¬ 
lations of fd/h{t) are not present when the semi-major axis is 
used to measure the disk fraction instead the instantaneous 
distance of the particles from the star. 

In Fig. we show the dependence of Tkl on the peri- 
center of the planetary orbit Rp (i.e., on es and Pb) for 
different semi-major axes of the planetesimals Pd (i-e., dif¬ 
ferent Pd). Tkl for Rp between 1 AU and 20 AU ranges from 
about 0.004 to 1 Myr depending on Pd. The wobbles happen 
on the timescale of ~ 0.1 Myr which is generally consistent 
(considering the factor a) with the Tkl for the particles in 
the inner parts of the disk and pericenter Rp ~ 1-5 AU and 
the full radial range of the disk for larger Pp. 

We suggest that the combination of the perturbation 
of the planetesimals orbits and a mechanism similar to the 
Kozai-Lidov oscillations leads to wobbling of the disk, when 
the eccentricities, inclinations, and the argument of periapse 
(i.e., the orientation of the orbits) change for a number of 
disk planetesimals. We illustrate the process in Fig. where 
we show snapshots of the simulation with the planetary peri¬ 
center at lAU and the disk inclination of 45°. The four 
snapshots show the initial state of the system, the times 
close to the minima {t = 0.3 and 0.9 Myr) and maximum 
{t = 0.6Myr) of the fd/h{t) modulation (see Fig. |^. At 
t = 0.3 and 0.9 Myr, the particles are collectively perturbed 
to higher inclinations and eccentricities and the plane of the 
disk is close to perpendicular to the orbital plane of the 
planet, while at t = 0.6 Myr, the disk has similar conhg- 
uration as in the beginning but with retro-grade rotation 
(inclination of about —45°). 



Figure 9. Timescale of the Kozai-Lidov mechanism, Tkl as 
given by Eq. as a function of the pericenter of the planetary 
orbit, Pp. Different lines show the dependence for different semi¬ 
major axes of the disk planetesimals, Pj. Several values of Pj 
are indicated in the plot. The dashed lines show the cases when 
Pb is outside the initial disk, while the full lines show the cases 
within the initial disk with a step-size of 20 AU. 


4.2 Short-term oscillations of fd/h 

Apart from the modulation on the timescales of ~ 0.1 Myr, 
the disk fractions fd/b{t) show periodical modulation with 
amplitudes of 0.03 and timescales of 0.05 Myr for most 
of the configurations (see Figs. andespecially the cases 
with higher disk fractions). The modulation results from res¬ 
onant spiral density waves and rings induced by the planet 
in the disk. If a resonant radius is located close to the initial 
outer edge of the disk, certain number of planetesimals orbit 
periodically just inside or outside the disk. The modulation 
is most prominent for the cases when the relative mass is 
fd/h ~0.7 and the resonant patterns are stable enough. If 
such resonant features are resolved by future observations, 
they can provide constraints on the orbit of the planet. 
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5 CONCLUSIONS 


We studied the lifetime of the debris disk in the peculiar 
system HD 106906. This ISMyr old star is orbited by a de¬ 
bris disk and a planetary mass companion at a separation of 
650 AU. We carried out simulations of the system using the 
AMUSE environment. Since the disk is much less massive 
than the star or the planet, we represent its planetesimals 
by zero-mass particles. We implemented a hybrid numerical 
method in which the orbit of the planet is solved indepen¬ 
dently of the disk and the disk planetesimals are integrated 
in the potential of the star and the planet. The initial condi¬ 
tions for the simulations were given by the observed charac¬ 
teristics of the system and the unconstrained characteristics 
of the system — namely the pericenter distance of the plan¬ 
etary orbit and the inclination of the disk with respect to 
the planetary orbit — were systematically varied. 

We find that more than 80% of the disk particles stay 
bound to the star for majority of the considered configu¬ 
rations and only in the case of orbits with low inclination 
10° and pericenter of the planetary orbit ^ 50 AU, a sub¬ 
stantial part of the disk is lost during the first 1 Myr of the 
evolution. To estimate how long the disk stays in a con¬ 
figuration consistent with the observations, we followed the 
ratio of the number of the disk particles with distance within 
the constrained disk radii (20-120 AU) and the number of 
the particles bound to the system. We define the lifetime 
of the disk when more particles are orbiting outside than 
within the constrained disk radius (i.e., more particles have 
is at distance < 20 AU or > 120 AU from the star). The 
lifetime of the disk is shorter than 1 Myr for orbits with low 
inclination i < 5° and comparable with 1 Myr when i ~ 5- 
10°, and with pericenter smaller or close to the inner edge 
of the disk {Rp ,$,50AU, see Figs. and [^. Such orbits 
are expected in the case when the planet formed closer to 
the star, most probably within the inner disk edge where 
it cleared the inner region, and was scattered to its cur¬ 
rent orbit by other member of the system. However, such 
interaction is estimated to occur during the first 10 Myr of 


the lifetime of planetary systems (e.g., Veras et al. 2009 
Scharf fc Menou||2009[|. Considering the current age of the 
system of 13 ± 2 Myr ( [Pecaut et al.|2012| ), we conclude that 
the configurations with lifetimes shorter than 1 Myr (i ^ 10° 
and Rp .$50AU) are inconsistent with the scenario accord¬ 
ing to which the current orbit resulted from planet-planet 
scattering from the inner disk. 

When the disk is inclined with respect to the plane¬ 
tary orbit with inclination ^1,40°, it can survive longer than 
1 Myr even in case the pericenter is within the inner disk 
edge. In these configurations, the disk wobbles (see Fig. |^. 
We argue that this is caused by a mechanism similar to the 
Kozai-Lidov oscillations induced by the planet on the disk 
particles. The planet can also induce resonant structures in 
the disk, such as spiral arms and rings. 


ACKNOWLEDGMENTS 

We thank the anonymous referee for reviewing our work 
and for insightful comments which improved the manuscript. 
We thank Guilherme Gongalves Ferrari, Inti Pelupessy, and 
Tjarda Boekholt for their help with the Kepler solver used 


in our numerical method. We thank Adrian Hamers for 
enriching discussions about the Kozai-Lidov mechanism, 
and Matthew Kenworthy and Tiffany Meshkat for discus¬ 
sions about HD 106906. We thank Grainne Gostigan for 
reading the manuscript and for her useful comments. This 
work was supported by the Interuniversity Attraction Poles 
Programme initiated by the Belgian Science Policy Of¬ 
fice (lAP P7/08 GHARM) and by the Netherlands Re¬ 
search Council NWO (grants #643.200.503, #639.073.803 
and #614.061.608) and by the Netherlands Research School 
for Astronomy (NOVA). 


REFERENCES 

Andrews S. M., Williams J. P., 2007, ApJ, 659, 705 
Bailey V., Meshkat T., Reiter M., Morzinski K., Males J., 
Su K. Y. L., Hinz P. M., Kenworthy M., Stark D., Mama- 
jek E., Briguglio R., Close L. M., Follette K. B., Puglisi A., 
Rodigas T., Weinberger A. J., Xompero M., 2014, ApJ, 
780, L4 

Baruteau C., Crida A., Paardekooper S.-J., Masset F., 
Guilet J., Bitsch B., Nelson R. P., Kley W., Papaloizou J. 
C. B., 2013, ArXiv e-prints, 1312, 4293 
Boley A. C., Payne M. J., Ford E. B., 2012, ApJ, 754, 57 
Boss A. R, 2011, ApJ, 731, 74 

Chabrier G., Johansen A., Janson M., Rafikov R., 2014, 
ArXiv e-prints, 1401, 7559 

Chatterjee S., Ford E. B., Matsumura S., Rasio F. A., 2008, 
ApJ, 686, 580 

Chen C. H., Jura M., Gordon K. D., Blaylock M., 2005, 
ApJ, 623, 493 

Chen C. H., Mamajek E. E., Bitner M. A., Pecaut M., Su 
K. Y. L., Weinberger A. J., 2011, ApJ, 738, 122 
Chiang E., Kite E., Kalas P., Graham J. R., Clampin M., 
2009, ApJ, 693, 734 

Crida A., Masset F., Morbidelli A., 2009, ApJ, 705, L148 
D’Angelo G., Durisen R. H., Lissauer J. J., 2011, in , Exo¬ 
planets, edited by S. Seager. Tucson, AZ: University of 
Arizona Press, 2011, 526 pp. ISBN 978-0-8165-2945-2., 
p.319-346. pp 319-346 

Davies M. B., Adams F. C., Armitage P., Chambers J., 
Ford E., Morbidelli A., Raymond S. N., Veras D., 2013, 
ArXiv e-prints, 1311, 6816 

de Zeeuw P. T., Hoogerwerf R., de Bruijne J. H. J., Brown 
A. G. A., Blaauw A., 1999, AJ, 117, 354 
Faramaz V., Beust H., Thebault P., Augereau J.-C., Bon- 
sor A., del Burgo C., Ertel S., Marshall J. P., Milli J., 
Montesinos B., Mora A., Bryden G., Danchi W., Eiroa 
C., White G. J., Wolf S., 2013, ArXiv e-prints, 1312, 5146 
Fujii M., Iwasawa M., Funato Y., Makino J., 2007, PASJ, 
59, 1095 

GonQalves Ferrari G., Boekholt T., Portegies Zwart S. F., 
2014, MNRAS, 440, 719 

Hamers A. S., Pols O. R., Claeys J. S. W., Nelemans G., 
2013, MNRAS, 430, 2262 

Hao W., Kouwenhoven M. B. N., Spurzem R., 2013, MN¬ 
RAS, 433, 867 

Hayashi C., 1981, Progress of Theoretical Physics Supple¬ 
ment, 70, 35 

Innanen K. A., Zheng J. Q., Mikkola S., Valtonen M. J., 
1997, AJ, 113, 1915 









10 Lucie Jilkovd and Simon Portegies Zwart 


Kozai Y., 1962, AJ, 67, 591 

Kraus A. L., Ireland M. J., Cieza L. A., Hinkley S., Dupuy 
T. J., Bowler B. R, Liu M. C., 2014, ApJ, 781, 20 
Lafreniere D., Jayawardhana R., van Kerkwijk M. H., 2008, 
ApJ, 689, L153 

Larwood J. D., Kalas P. G., 2001, MNRAS, 323, 402 
Lestrade J.-F., Morey E., Lassus A., Phou N., 2011, A&A, 
532, 120 

Lidov M. L., 1962, Planet. Space Sci., 9, 719 
Malmberg D., Davies M. B., Heggie D. C., 2011, MNRAS, 
411, 859 

Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398 
Naud M.-E., Artigau E., Malo L., Albert L., Doyon R., 
Lafreniare D., Gagne J., Saumon D., Morley G. V., Allard 
F., Homeier D., Beichman C. A., Gelino C. R., Boucher 
A., 2014, ApJ, 787, 5 

Pearce T. D., Wyatt M. C., 2014, ArXiv e-prints 
Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, ApJ, 746, 
154 

Pelupessy F. L, Janes J., Portegies Zwart S., 2012, New As- 
tron., 17, 711 

Pelupessy F. L, van Elteren A., de Vries N., McMillan S. 

L. W., Drost N., Portegies Zwart S. F., 2013, A&A, 557, 
84 

Perets H. B., Kouwenhoven M. B. N., 2012, ApJ, 750, 83 
Portegies Zwart S., Boekholt T., 2014, ApJ, 785, L3 
Portegies Zwart S., McMillan S., Harfst S., Groen D., Fuji! 

M. , Nuallain B. O., Glebbeek E., Heggie D., Lombardi J., 
Hut P., Angelou V., Banerjee S., Belkus H., Fragos T., 
Fregeau J., Gaburov E., Izzard R., Juric M., Justham S., 
Sottoriva A., Teuben P., van Bever J., Yaron O., Zemp 
M., 2009, New Astron., 14, 369 

Portegies Zwart S. F., McMillan S. L. W., 2005, ApJ, 633, 
L141 

Ralikov R. R., 2011, ApJ, 727, 86 

Raymond S. N., Armitage P. J., Gorelick N., 2010, ApJ, 
711, 772 

Rebull L. M., Stauffer J. R., Megeath S. T., Hora J. L., 
Hartmann L., 2006, ApJ, 646, 297 
Reche R., Beust H., Augereau J.-C., Absil O., 2008, A&A, 
480, 551 

Scharf G., Menou K., 2009, ApJ, 693, L113 
Steinhausen M., Olczak C., Pfalzner S., 2012, A&A, 538, 
AlO 

Thebault R, 2012, A&A, 537, 65 

van Leeuwen F., 2007, A&A, 474, 653 

Veras D., Crepp J. R., Ford E. B., 2009, ApJ, 696, 1600 

Veras D., Evans N. W., 2013, MNRAS, 430, 403 

Vorobyov E. L, 2013, A&A, 552, 129 

Wisdom J., Holman M., 1991, AJ, 102, 1528 

Wyatt M. C., 2003, ApJ, 598, 1321 

Wyatt M. C., 2008, Annual Review of Astronomy and As¬ 
trophysics, 46, 339 



